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ABSTRACT 

P-] ■ We present an analysis of the high frequency X-ray variability of NGC 4051 (Mbh ~ 1.7 x 

10 6 Mq) based on a series of XMM-Newton observations taken in 2009 with a total exposure 
of ~ 570 ks (EPIC pn). These data reveal the form of the power spectrum over frequencies 
from 10~ 4 Hz, below the previously detected power spectral break, to ;> 10 -2 Hz, above 
the frequency of the innermost stable circular orbit (ISCO) around the black hole (Visco ~ 
10~ 3 — 10 -2 Hz, depending on the black hole spin parameter j). This is equivalent to probing 
' frequencies J> 1 kHz in a stellar mass (Mbh ~ 10 Mq) black hole binary system. The 

power spectrum is a featureless power law over the region of the expected ISCO frequency, 
suggesting no strong enhancement or change in the variability at the fastest orbital period in 
the system. Despite the huge amplitude of the flux variations between the observations (peak- 
to-peak factor of ;> 50) the power spectrum appears to be stationary in shape, and varies 
in amplitude at all observed frequencies following the previously established linear rms-flux 
! relation. The rms-flux relation is offset in flux by a small and energy-dependent amount. 

The simplest interpretation of the offset is in terms of a very soft spectral component that 
is practically constant (compared to the primary source of variability). One possible origin 
for this emission is a circumnuclear shock energised by a radiatively driven outflow from the 
central regions, and emitting via inverse-Compton scattering of the central engine's optical- 
UV continuum (as inferred from a separate analysis of the energy spectrum). A comparison 
with the power spectrum of a long XMM-Newton observation taken in 2001 gives only weak 
evidence for non-stationarity in power spectral shape or amplitude. Despite being among the 
most precisely estimated power spectra for any active galaxy, we find no strong evidence for 
quasi-periodic oscillations (QPOs) and determine an upper limit on the strength of a plausible 
/\ ' QPO of <; 2 per cent rms in the 3 x 10~ 3 — 0.1 Hz range, and ^5 — 10 per cent in the 

^ . 10~ 4 — 3 x 10~ 3 Hz range. We compare these results to the known properties of accreting 

stellar mass black holes in X-ray binaries, with the aim further of developing a 'black hole 
unification' scheme. 

Key words: galaxies: active - galaxies: individual: NGC 4051 - galaxies: Seyfert - X-rays: 
galaxies 



1 INTRODUCTION 

If strong gravity dominates the dynamics of the inner accre- 
tion flows around black holes then an elementary consequence 
is scale invariance: many important aspects of accretion onto su- 
permassive black holes (Mbh J> 10 6 Mq) in Active Galac- 
tic Nuclei (AGN) should be fundamentally the same as for stel- 
lar mass black holes (Mbh ~ 10 Mq) in Galactic B lack 
Hole systems (GBHs) . (See e.g. IShakura & Sunvaevl 1 19761 and 
Mushotzky et al. I ll993l .) Over the past few years several similari- 
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ties have been observed in the X-ray variability of near by AGN and 
GBHs, supporting the idea of black hole unification dFender et al.l 
2006). In particular, the power spectrum (often power spectral 
density, PSD), and inter-band X-ray time lags appear to have 
similar frequency dependence and relative amplitudes in AGN 
and GBHs but with characteristic frequencies sc aled (to first or- 
der) inversely with mass, v oc 1/Mbh (see e.g. lLawrence et al 



1987t lHavashida et alJl998tlEde lson & Na ndrall 19991; lUttlev et al 



l2002l ; lyaughan et alj|20 03: Markowitz et al. 2003; M c Har dy et al 
| 2004l; iDone & Gierlinskil 120051: hyTHardv et al.l I2006L 120071 : 
Gierliriski e t al.l2008l : lArevalo et alj2008l) . Another property of the 
X-ray variability common to both types of system is a linear rms- 
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Figure 1. Time series of 0.2 — 10 keV EPIC pn count rate (background-subtracted) from each of the 15 observations taken in 2009. The bin size is At = 100 
s. The 15 separate ~ 30 — 40 ks observations, taken typically 2 — 4 days apart, have been concatenated in this plot, to allow easy comparison between them, 
with the end points marked by the dotted vertical line. The background count rate is also shown as a grey curve, but is almost always far weaker than the 
source. The labels above each section of data refer to the spacecraft revolutions during which the observations were made. 



flux relation, a scaling of the rms amplitude of variability with av- 
erage sourc e flux that appears to hold over a very wide range of 
timescales dUttlev & M c Hardvll200ll; Vaug han et al]|2003l ;lGaskell 
l2004LlGleissner et al.ll2004l:IUttlev et al.ll2005h . 

Arguably among the most important recent breakthroughs 
in the study of stellar mass black holes was the discov- 
ery of quasi-periodic oscillations (QP Os) at high frequen- 
cies (v Q > 1 00 Hz) in GBHs (see iRemillard eTai] 1 19991 ; 



Strohmavenl2001l: IRemillard et alJl2002l. 120031 : Ivan der Klisll2006l : 
Remillard & McClintocla2006t) . These are the fastest variations ob- 
served from GBHs and their high frequencies, close to the Ke- 
plerian frequency of the innermost stable circular orbit (ISCO; 
lAsco), indicate an origin in the central regions of the accretion 
flow. High frequency QPOs (HF QPOs) therefore carry informa- 
tion on the strongly curved spacetime close to the event horizon and 
should, in principle, constrain the fundamental parameters of black 
holes: mass and spin. As a product of strong gravity, HF QPOs 
should also be present in AGN, but until recently there were no ro- 
bust detections, due mainly to t he insufficient length of previou s 
XMM-Newton observations (see IVaughan & Uttlevl 1200a . l2006h - 
Recently a QPO was d etected in the Seyfert galaxy RE J 1034+396 
dGierlinski et alj| 2008 ) . This is a major discovery that adds further 
support to the general idea that AGN are scaled-up GBHs, but this 
one detection provides limited diagnostic power because it is so far 
unique and the other proper t ies of the source are not well under- 
stood dMiddleton et al.ll2009l : lMiddleton & DoTiehoiOh 

In parallel with these X-ray advances, studies of transient ra- 
dio jets from GBHs have led to the development of a scheme 
that unifies jet production with accretion state dFender et al.l 
12004b . Relativistic jets must be launched from the region dom- 
inated by strong-field gravity, and so scale invariance implies 
this accretion state/jet unification scheme should extend to AGN 



Heinz & Sunyaev 2003), a hypothesis supported by observations 



(Merl oni etalj20 03; Falcke et al. 2004; Kor ding et al.l2006j) . From 



tion for both GBH "states" and AGN "types" may be unified into a 
single framework for the activity cycles of accreting black holes. 

Here we discuss the results obtained from a series of XMM- 
Newton observations of the bright, highly- variable Seyfert 1 galaxy 
NGC 4051. These observations were designed to provide some 
of the best constraints to date on the high-frequency X-ray vari- 
ability properties of any AGN and improve our understanding of 
black hole unification generally. This particular target was cho- 
sen because of a combination of fa ctors. NGC 40 5 1 (redshift z = 
0.002336, distance D « 15 Mpc: lRusselill2002l) has a relatively 
low and well-de termined black hole mass, A/bh ~ 1-7 ± 0.5 x 
10 6 Mq dDennev et al]|2009l I20I0I) . disp lays strong, persistent 
X-ray flux and spectral variability (see e.g. iLawrence et al.|[l987l: 
Papad akis & Lawrence] Il995l: llVTHardv et al.| 120041 : lUttlev et al.l 
120041 : IPonti et alJl2006l ; TTerashima et al.ll2009h . an d a soft X-ray 



spect r um rich in em i ssion a nd absorption feat ures (Collinge et al 



these results a new paradigm is emerging in which accretion mode, 
X-ray spectrum, high-frequency timing properties and jet produc- 



| 200ll : IPounds et all |2004| ; bgle et alj |2004| ; Isteenbrugge etal 
1200' 



2 OBSERVATIONS AND DATA ANALYSIS 

NGC 4051 was observed by XMM-Newton 15 times over 45 days 
during May- June 2009, and previously during 2001 and 2002 (see 
table [T}. The total duration of the useful data from the 2009 cam- 
paign is ~ 572 ks, giving ~ 6 x 10 6 EPIC pn source counts (us- 
ing PATTERN — 4) or ~ 9 x 10 6 EPIC source counts (using 
PATTERN - 12 pn and MOS combined). 

The raw data were processed from Observation Data Files 
(ODFs) following standard procedures using the XMM-Newton 
Science Analysis System (SAS vlO.0.2). The EPIC data were pro- 
cessed using the standard SAS processing chains to produce cal- 
ibrated event lists. For each observation, source events were ex- 
tracted from these using a 40 arcsec circular region centred on the 
the target, and background events were extracted from a rectangular 
region on the same chip but not overlapping with the source region. 
(In order to obtain a precise background estimate the background 
region was ~ 7 times larger than the source region.) Examination 
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Table 1. Observation log. The columns list (1) the observation identifier, (2) 
the spacecraft revolution number, (3) the date of the start of the observation, 
(4) the observation duration. 



Obs ID. 




Obs. date 


Duration 


ID 


rev. 


(YYYY-MM-DD) 


(ks) 


0109141401 


0263 


2001-05-16 


121958 


0157560101 


0541 


2002-11-22 


51866 


0606320101 


1721 


2009-05-03 


45717 


0606320201 


1722 


2009-05-05 


45645 


0606320301 


1724 


2009-05-09 


45548 


0606320401 


1725 


2009-05-11 


45447 


0606321301 


1727 


2009-05-15 


32644 


0606321401 


1728 


2009-05-17 


42433 


0606321501 


1729 


2009-05-19 


41813 


0606321601 


1730 


2009-05-21 


41936 


0606321701 


1733 


2009-05-27 


44919 


0606321801 


1734 


2009-05-29 


43726 


0606321901 


1736 


2009-06-02 


44946 


0606322001 


1737 


2009-06-04 


39756 


0606322101 


1739 


2009-06-08 


43545 


0606322201 


1740 


2009-06-10 


44453 


0606322301 


1743 


2009-06-16 


42717 



0.1000 - 
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Figure 2. Power spectrum estimate (histogram) for the broad band (0.2—10 
keV) EPIC pn data, computed using 10 ks intervals at Ai = 5 bins, with 
best-fitting bending power law continuum model (solid curve). Also shown 
(grey circles) are two low frequency power density estimates computed us- 
ing 28 ks intervals. These lower frequency data were not used in the fitting 
(see text for details). 



of the background time series showed the background to be rela- 
tively low and stable throughout most of the observations, with the 
exception of the final few ks of each observation where the back- 
ground rate increased as the spacecraft approached the radiation 
belts at perigee (each of the observations occurred towards the end 
of a spacecraft revolution). The EPIC data were taken using the 
small window mode of each camera and the medium blocking fil- 
ter to reduce pile-up and optical loading effects, respectively. Nev- 
ertheless, during bright intervals NGC 4051 is sufficiently bright 
to cause moderate pile-up, espe cially in the MPS, w hich can af- 
fect power spectrum estimation iTomsick et alj [2004). In order to 
mitigate any distorting effects all results are based on pn data PAT- 
TERNS — 4 only, but where possible checked for consistency with 
the MOS data. Response matrices were generated using RMFGEN 
vl . 55 . 2 and ancillary responses were generated with ARFGEN 
v. 1 . 77 . 4. 

Time series were extracted from the filtered events for the 
source and background regions using a range of time bin sizes and 
energy ranges. These were all corrected for exposure losses, as cat- 
alogued in the Good Time Information (GTI) header contained in 
each event file, except where less than 30 per cent of a bin was 
exposed, in which case the data were calculated by linear inter- 
polation from the nearest 'good' data either side and appropriate 
Poisson noise was added (this was needed for only 0.3 per cent 
of the total exposure time from the pn, for At = 5 s time bins). 
The outputs of this process were regularly sampled and uninter- 
rupted, exposure-corrected time series for source and background. 
The background time series was then subtracted from the corre- 
sponding source time series (after scaling by the ratio of the 'good' 
extraction region areas). Figure [T] shows the resulting time series 
from the 15 separate 2009 observations. 



3 POWER SPECTRUM 

The power sp ectrum of the 2009 data was estimated using stan- 
dard methods ( van der Kfiiri989). In particular, time series were 
extracted with binning At = 5 s and divided into uninterrupted 



10 ks intervals. For eac h interval a periodogram was computed in 
units of relative pow er ( Mivamot oeTafll 199ll ; Ivan der Klislll997l : 
Vaughan et al. 2003), the Poisson noise leveQ was subtracted, and 
then the periodgrams were averagec0. Figure|2]shows the resulting 
power spectrum estimate after rebinning in logarithmic frequency 
intervals. The model fitting was performed using XSPEC v 12.6.0 
dAmaud[l996l) using x 2 as the fit statistic, all data were binned to 
ensure that Gauss-normal statistics applied, and errors quoted cor- 
respond to 90 per cent confidence limits (e.g. using a A% 2 = 2.71 
criterion for one parameter) unless otherwise stated. 

The previous ana lysis of the power spectrum of NGC 405 1 by 
|M C Hardy et al.l |2004) found a smoothly bending power law con- 
tinuum provided an adequate description of the data 
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with indices of aa ow ~ —1.1 at low frequencies and ahigh ~ — 2 
above fbcnd ~ 8 x 10 -4 Hz. This model provided an excel- 
lent fit to the 2009 power spectrum, as shown in Figure [2] with 
the lower frequency index fixed at ai ow = —1.1 (since it i s 
much better constrained by the RXTE data; lM c Hardv et al.ll2004T) . 
The best-fitting values for the free parameters of the model were 
Qhi g h = -2.31 ± 0.08 and ^ bcn d = 2.3 ± 1.0 x 10~ 4 Hz. The fit 
statistic was % 2 = 31.19 with 27 degrees of freedom (dof), giving 



1 The expected power density due to measurement uncertainties in this case 
is Pjv = 2AT( cr 2 ), where AT is the sampling rate and (a 2 ) is the mean 
square error. See lVaughan et alj fc003l) , 

2 We have confirmed that the power spectral estimation and modelling 
was not significantly biased by sampling distort ions - commonly known 
as 'ali a sing' and 'red noise l eakage,' see lvan der Klij fl989h ; luttlev et alj 
1 2002); Vaughan et al. (2003). Aliasing is minimal for contiguously binned 
data Ivan der Klis 1989), like those used here, and does not occur for the 
Poisson noise that dominates the highest frequencies probed. Leakage could 
in principle have been a problem, resulting from the finite duration of each 
observation, but in practice was not significant because the power spectrum 
was sufficiently flat (~ v~ ) on timescales comparable to the observation 
length. This was confirmed using an analysis discussed in section [4] foot- 
note!?] 
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a p- value of 0.26. Allowing ai ow to vary provided no improvement 
in the fit (Ax 2 < 2) and the index itself was poorly constrained 
(with a 90 per cent confidence interval of [—0.25, —1.78]), and 
restricting «i ow to vary only wit hin the range [—1.2,- 0.91 (the 
confidence interval determined by M c Hard v et alj J2004) from the 
RXTE data) provided very little change in the confidence interval 
for Vhend compared to the fixed ai ow = —1.1 model. There is 
no obvious indication, in either the overall fit statistic or the resid- 
ual plot, of additional power spectral components. We discuss this 
point in more detail in Section[5] 

A sharply broken power law model also provided a good fit, 
with x 2 = 37.52 for 27 dof (p = 0.086), slightly worse than 
the smoothly bending model. The index below the break was again 
fixed to aiow = —1.1, and the best-fitting break frequency was 
^Wak = 2.5 ± 0.6 x 10~ 4 Hz, and a high frequency index of 
Ohigh = —2.18 ± 0.05. These parameters, and the overall fit qual- 
ity, are very similar to those of the bending power law model, indi- 
cating that the data are rather insensitive to the detailed shape of the 
change from flat to steep index. By contrast, the power spectrum is 
less well explained in terms of an exponentially cut-off power law, 
which gives x 2 = 44.03 for 26 dof (p = 0.015), but only by allow- 
ing Qio w ~ —1.8 which is i nconsistent with the long-term RXTE 
results dM c Hardv et afll20Q4 Assuming qi ow = 1.1 gave an un- 
acceptable fit (x 2 = 217.81 for 27 dof; p < 0.001). 

The model fitting was repeated using a Bayesian scheme by 
treating the likelihood function as ~ exp(— x' 2 / 2) and assignin g 
simple prior distributions to the parameters. See Vaughan (2010), 
and references therein, for a brief introduction to these ideas. Uni- 
form priors were assigned to a^igh and the overall normalisation, 
and to log ^bcnd (equivalent to a p(f) = 1/v Jeffreys' prior), al- 
though in practice using uniform or Jeffreys' prior made no sig- 
nificant difference to the results. As expected, the parameter val- 
ues at the posterior mode were almost identical to those giving the 
best-fit using the classical minx 2 method. But having the prob- 
lem restated in Bayesian terms allowed for an exploration of the 
posterior distribution of the parameters, p(8c | x obs ) with 9c the 
parameters specifying the continuum and x obs the observed data, 
using a Markov Chain Monte Carlo (MCMC) scheme to randomly 
draw sets of parameters from the posteriori The marginal density 
for the bend frequency, pf^cnd | x obs ), as computed using the 
MCMC data, is shown in Figure[3j the 90 per cent credible interval 
on ^bcnd was [1.3, 3.2] x 10 -4 Hz. 

The limited durations of the individual 2009 observations 
makes it difficult to directly probe frequencies lower than 10 -4 
Hz. But for completeness the power spectrum was also estimated 
using longer uninterrupted intervals (28 ks duration, of which there 
were 15 in 2009, one for each observation). However, the estimates 
of power density at the lowest frequencies were not normally dis- 
tributed in this cas e, due to the small number of e stimates contribut- 
ing to the average JPapadakis & Lawrencdl 19931) . These data were 
not included in the power spectral fitting, but the lowest frequency 
points are shown in Figure [2] and are clearly consistent with an ex- 
trapolation of the best-fitting model. 

As a test for a possible high frequency cut-off in the power 



3 Specifically, five chains of length 5 X 10 4 were generated (after removing 
an equal length "burn-in" period), using a Metropolis-Hastings algorithm 
with a multivariate Cauchy proposal distrib ution. Convergence of the chains 
was confirmed using the R statistic i Gelman & Rubin 1992; Gelman et al. 
|2004|) . The results from each chain were then combined to produce a large 
sample of posterior draws. 
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Figure 3. Marginal posterior probability density for the bend frequency in 
the bending power law model, computed by MCMC (see text for details). 
That for the 2009 data is shown as a solid curve, and that for the 2001 data 
is shown as the dashed curve. 



spectrum the bending power law model was modified by includ- 
ing an exponential cut-off (the highecut model in XSPEC). This 
provided only a very small improvement in the overall fit statis- 
tic (A^ 2 = 4.1 for 2 additional free parameters) with only poorly 
constrained parameters (cut-off frequency i/ C ut ^ 8 x 10" 3 and 
e-folding frequency Vf \d = 2 ± 1 x 10 -2 ). The power spectrum 
estimated from the high flux data (see next section) has slightly 
higher signal-to-noise above 10 mHz, but again shows no signifi- 
cant difference in the quality of the fit between the bending power 
law with and without an exponential cut-off at high frequencies. 



3.1 Energy dependence of the power spectrum 

The dependence of the power spectrum on photon energy was in- 
vestigated by dividing the data into three broad energy bands - 
0.2-0.7, 0.7-2 and 2-10 - and estimating and fitting the power 
spectrum for each of these. A simultaneous fit to all these spectra 
using the same continuum model (equation QJ gave a rather poor 
fit (x 2 = 118.9 with 85 dof, p = 0.009) with the parameters tied 
between the three spectra. But allowing to be different for each en- 
ergy band, while keeping I4, e nd tied between them, gave a much 
better fit (x = 98.5 with 83 dof, p — 0.12); the improvement is 
significant in an F-test, with p < 0.001, indicating the high fre- 
quency power spectrum is energy dependent. The slopes estimated 
for the low, medium and high energy bands were 2.31 ± 0.07, 
2.24 ± 0.09 and 2.04 ± 0.10, respectively. The flattening of the 
power spectrum at hig her energies is very similar to that detected by 
lM c Hardv et al.l (T2004) from the 2001 d ata, and previously observed 
in several other Seyfert galaxies (e.g. iNand ra & Papadakis l200ll ; 
I Vaughan & FabiarJ 120031) . Allowing the bend frequency Vbend to 
be different for each energy band, but with a n i g h tied between 
them, provided a considerably worse fit (x 2 = 108.0 with 83 
dof) than the energy-dependent ahigh model. Furthermore, allow- 
ing both i^bend and «hi gn to differ between energy bands did not 
improve significantly on the variable augh fit (x 2 = 95.7 with 
81 dof). In this case the best-fitting fbead values for the low and 
medium energy bands were very similar while the hard band gave 
a somewhat lower value, but with a much larger confidence inter- 
val. We conclude there is strong evidence for a change in the high 
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Figure 4. Power spectrum estimates computed as in Figure [2] from 10 ks 
intervals but averaged into three flux bins, with average (0.2 — 10 keV 
pn) count rates of 4.8, 8.4 and 17.4 ct s -1 . The lower panel shows the 
ratio of highest to lowest flux power spectral estimates, which is consistent 
with a constant over two decades in frequency, indicating a strong increase 
in absolute power but no strong change in shape with average flux. The 
data were binned logarithmically in frequency, by a factor 1.4, and the bin 
centres shifted slightly to prevent the error bars overlapping. 



frequency slope with energy, but the data do not resolve any differ- 
ence in the location of the bend. 



4 RMS-FLUX RELATION 

The 15 observations shown in Figure[T]clearly show very different 
mean levels and amplitudes of variability, with noticeably stronger 
variations occurring preferentially during brighter intervals (e.g. 
revl730), and the lowest flux intervals being extremely stable (e.g. 
rev 1739). This behaviour i s exactly as expected for sources that 
obey the rms-flux relation dUttlev & M c Hardvll200ll ; lUttlev et all 
l2005t) . 

In order to test this explicitly and in detail the power spectrum 
was estimated in different flux bins. As before, the time series were 
divided into uninterrupted 10 ks intervals, the periodograms com- 
puted and the contribution of the Poisson noise subtracted, but us- 
ing absolute u nits of power (i.e. w ithout normalising to the squared 
mean rate, see lVaughan et alj|2003h . The intervals were sorted ac- 
cording to their mean flux, and the periodograms were then aver- 
aged in flux bins to provide estimates of the power spectrum as a 
function of mean flux. Figure[4]shows the power spectral estimates 
when the data were divided into three flux bins. There is clearly 
a large difference in the overall amplitude (in absolute units) but 
little change in shape, as revealed by the ratio of bright and faint 
power spectra. Direct fitting of the high and low flux power spectra 
gives the same result: the best-fitting values for Qhigh an d ^bend es- 
timated for the high flux data fall within the 90 per cent confidence 
intervals derived fr om the low flux dat a . The same effect was seen 
in Cygnus X-l by lUttlev & tvTHardvl J200lh (see their figure 2), 
and is consistent with there being a single rms-flux relation inde- 
pendent of the frequencies used to compute the rms. 

The detailed rms-flux relation was then examined explicitly 
by dividing the data into 1 ks uninterrupted intervals, computing 
a periodogram for each (again, in absolute units), subtracting the 
Poisson noise from each, and then averaging the results in several 



Figure 5. 0.2 — 10 keV rms-flux relation. The rms is measured over the 
1 — 10 mHz range (see text for details). The solid line marks the best-fitting 
linear model. The lower panel shows the standardised residuals, i.e. (data - 
model)/error. 

fluxbimQ (Fi). The 1 — 10 mHz rms, Oi was then estimated for each 
flux bin by integrating the corresponding power spectrum estimate 
over this frequency range. The resulting rms-flux data are shown in 
Figure [5] (The uncertainties on the rms estimates were calculated 
based on the variance of the sum of the averaged periodograms, 
which are themselves calculated based on the chi-square distribu- 
tion of periodogram ordinates. The method is explained more fully 
in Heil, Vaughan & Uttley in prep.) 

The dependence of rms on flux is clearly very close to lin- 
ear over the full flux range, as previousl y noted in this source 
and accreting blac k holes generally (e.g. lUttlev & M c Hardvll200ll; 
lUttlev etaf. 2005, Heil, Vaughan & Uttley in prep.). The best-fitting 
linear function - Oi = k(Fi — C), where k is the gradient and 
C is the offset on the flux axis - gave a rather high fit statistic 
(X 2 = 30.05 for 12 dof) but the residuals show little structure. 
The best-fitting linear function has a negative offset in the rms axis, 
equivalent to a positive offset on the flux axis. The energy spec- 
trum of this offset was computed by dividing the pn data into dif- 
ferent energy bands (14 logarithmically spaced) and for each one 
calculating the rms-flux data and fitting with a linear model. The 
spectrum of the flux offset C(E) is shown in Figure(6] along with 
some representative pn spectra from bright (revl730), intermediate 
(revl743) and faint (revl739) flux observations. The spectra were 
"fluxed" for display purposes, that is, they have been normalised 
to flux density units by dividing out the effective area curve (see 
Appendix |A] for details). The spectrum of the offsets from the rms- 
flux relation is clearly soft: it is significantly positive in the lower 
energy bands, decreasing until consistent with zero intercept above 
~ 2 keV. Most interestingly, it has a very similar shape and strength 
to the faint flux spectrum below 2 keV. 



4 We have checked that the power spectra estimated using sho rt time series 
interv als (e.g. 1 ks) are not significantly biased by leakage (see lUttlev et alj 
12001 and references therein) in two different ways. First we compared the 
average power spectrum estimated from the real data using 1 ks, 10 ks and 
28 ks intervals and confirming they are all consistent. Secondly, by simu- 
lating data (using the best-fitting bending power law model) with the same 
sampling pattern as the real data, and estimating the bias between the power 
spectral estimate and the model. 
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IPottschmidt et al.1 J2003h for details of this parameterisation. In the 
present case we have assumed two different values, namely Q = 5 
and Q — 10, whic h are fairly typical values fo r high frequency 
QPOs in GBHs (e.g. lRemillard et alj|2002ll2003l) . Using these two 
representative values rather than allowing Q to be a free parameter 
made the QPO search and upper limit calculations (discussed later) 
much more efficient. 

One of the simplest ways to select between these two compet- 
ing models is to assess the improvement in the \ 2 fit statistic upon 
the addition of the QPO, i.e. the difference between the minimum 
X 2 for Ho and Hi . For a given dataset a large reduction (improve- 
ment) in the fit statistic between Ho and Hi is taken to be evidence 
in favour of a QPO. For the 2009 data, and the full 0.2 - 10 keV 
energy band, the reduction in \ 2 statistic was only 6.85 (assuming 
Q = 10, and a very similar value assuming Q — 5), with best- 
fitting parameters vq — 3.8 ± 0.3 mHz and R — 1.2 ± 0.5 per 
cent. The modest improvement in the fit falls short of our thresh- 
old for a significant detection, as discussed below. The QPO search 
was repeated for the three energy bands discussed above, with very 
similar results. The 2-10 keV band gave the largest improvement 
in the fit upon the addition of a QPO (A X 2 = 10.94 assuming 
Q = 5), but this still falls short of a detection in the a — 0.01 test. 

The lack of detection may itself be an interesting result if the 
observation was sufficiently sensitive, and to assess this requires a 
calculation of the upper limit on the strength of a QPO. The first 
step is to define a precise detection procedure and then calibrate 
its statistical 'power' as a function of QPO strength. That would 
allow us to determine the weakest QPOs that would be detected 
with reasonable probability. 



Energy (keV) 



Figure 6. "Fluxed" EPIC pn spectra shown in E X F(E) units. Moving 
from top to bottom the spectra are: (1) from a bright flux interval (rev 1730), 
(2) from a typical flux interval (rev 1743) (3) from a faint interval (rev 1739), 
(4) the offset of the rms-flux relation C(E) (see text). Errors for the rms- 
flux offset spectrum are the errors on the linear regression parameter. 



5 A SEARCH FOR QUASI-PERIODIC OSCILLATIONS 

The power spectral fitting discussed above provided no clear ev- 
idence for narrow features, i.e. QPOs. Such features are usually 
described in terms of Lorentzian "lines" seen in the power spec- 
trum in addition to the continuum. In order to assess the evidence 
for a QPO more rigorously we defined two competing hypotheses: 
the null hypothesis Ho of a continuum power spectrum, and the 
alternative hypothesis Hi which includes an additional Lorentzian 
component. 

H : P(y) = cont(v,6o) (2) 



Hi : P{v) = cont{v,6 c ) + qpo{v,e Q ). (3) 

where the continuum model is defined above (equation [T) and the 
QPO will be described by a Lorentzian profile 

2R 2 Qv/ir 



qpo{v,Q Q ) = 



(4) 



The parameters of the continuum are 9c = {cthigh, ^bcnd, N} and 
the QPO parameters are 8q — {vq, Q, R}, where vq is the cen- 
troid frequency, Q is the width parameter (Q — vq/Av, where 
Ai/ is the Full Width at Half Maximum), and R is the ampli- 
tude parameter (equal to the total rms in the limit of high Q). See 



5.1 Detection procedure 

The search for QPOs in binned power spectral data is no dif- 
ferent from the search for "lines" in any other one dimen- 
sional spectrum, and suffe rs from the same sta ti stical and com- 
putational challenges. S ee iFreeman et al.l jl999t) , IProtassov et al.l 
d2002l) , |Park et alj 120081) for elaboration of the issues. 

The difference in the minimum chi-square statistic for each 
model was used as a test statistic. For a dataset x we compute 



T(x) = X 2 (x, J fJ )-X (x,#i) 



(5) 



The nuisance parameters - specifying the continuum and the 
QPO - are eliminated by minimisation, i.e. x 2 ( x i^o) = 
ming c [x 2 ( x , Ho, 9c)\ This test statistic is closely related to 
the Likelihood Ratio Test (LRT), as has been discussed by 
IProtassov et al. (2002). A large value of the test statistic (i.e. a large 
improvement in the fit upon adding a QPO to the model) is taken 
as evidence for a QPO. 

The critical value of T for detection, T cr i t , was calculated by 
a Monte Carlo method as discussed in Appendix |B1 In the present 
case we found T cri t = 12.28 for a = 0.01 (this is for Q = 10, 
using Q = 5 gave 13.05 which is practically the same). In other 
words, assuming the continuum model to be true (with parame- 
ters represented by the posterior distribution p(9c \ x obs ) derived 
above from the real data), and no QPO, a reduction in the \ 2 fit 
statistic of T ^ 12.28 will occur with probability a — 0.01 (and 
have a posterior predictive p- value ^ 0.01). 
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Figure 7. Upper limit on the strength of a QPO, assumed to have a 
Lorentzian spectrum, defined for an a = 0.01 test. The solid curve was 
calculated assuming a Q = 10 QPO and the dotted curve was calculated 
assuming Q = 5. In each case the three curves show (from bottom to top) 
the upper limits corresponding = 0.5, 0.9 and > 0.995 powers. The 
detection power at each point in the plane was computed using 200 Monte 
Carlo simulations of data with a QPO present and recording the fraction of 
simulations that gave a detection (using a A\ 2 ^ 12.28 detection criterion 
for Q = 10). See section[3]and Appendix |S]f or details. 



5.2 Upper limit procedure 

In order to define an upper limit on the strength of possible QPOs 
we use the definition of an upper limit in terms of the powe^f] /3 of 
a detection procedure, a s discussed very clearly in the recent paper 
bv lKashvap eTail l |2010l) . A QPO of a given frequency and strength 
(vq,R) has a probability of detection P(vq,R) using the above 
method (for fixed a). 

We used a Monte Carlo method, as discussed in Appendix IB] 
to estimate the detection probability P(vq , R) for a range of differ- 
ent QPO locations and strengths (vq and R values). The estimate 
is simply the fraction of simulations (for each vq and R) that gave 
a significant detection defined by T(x slm ) ^ T cr it- Contours of 
constant /? on the vq-R plane define the upper limit on the strength 
J? of a QPO as a function of frequency vq. Figure|7]shows the con- 
tours for p = 0.5, 0.9 and > 0.995 assuming Q = 10 (solid 
curves) and Q — 5 (dotted curves). The middle of each set of 
curves represents the weakest QPO that could be present (at each 
frequency) and have a probability ^ 0.9 of being detected. This 
limit extends below R — 5 per cent for v > 3 x 10 -4 Hz, and 
below R = 2 per cent for v 3 x 10~ 3 Hz. The shape and ampli- 
tude of this curve is similar to that obtained usin g the simpler A\ 2 
confid ence contour method used previously by Vaughan & Uttlevl 
d2005b . 



6 LONG TERM STATIONARITY 

The power spectral estimation and fitting described in section [3] 
was repeated for the 2001 data. Figure [8] shows a comparison of 
the 2001 and 2009 power spectra which are in broad agreement 



Figure 8. Comparison of 2009 and 2001 power spectrum estimates. Top 
panel: the estimated power spectra computed for the 2001 observation (blue 
histogram) and the combined 2009 observations (black; same as Figure|2j. 
Bottom panel: ratio of the two estimates. 



in shape and overall (fractional) amplitude. The best-fitting bend- 
ing power law model (% 2 = 26.22 for 27 dof, p = 0.51) had a 
slightly higher bend frequency (u w 5.5 x 10~ 4 Hz with a 90 
per cent confidence interval of [3.0, 9.0] x 10 -4 Hz, and a cred- 
ible interval of [2.7,8.4] x 10~ Hz) and high-frequency index 
(ahigh = —2.48 ± 0.16) compared to the 2009 data. The marginal 
distribution of ^b cn d from both datasets are shown in Fig[3] These 
appear to suggest there was a slight change in the power spectral 
shape between the two sets of observations. It is for this reason that 
we have not combined the 2001 and 2009 data to produce a com- 
posite power spectrum. 



7 DISCUSSION AND CONCLUSIONS 

We have presented an analysis of the rapid X-ray flux variability 
of the low-mass Seyfert 1 galaxy NGC 4051 covering the range 
10~ 4 to ~ 2 x 10 -2 Hz, based on a series of 15 XMM-Newton ob- 
servations taken during 2009. This is among the best X-ray power 
spectra for any active galaxy, and the large range of fluxes sampled 
during the observations has allowed for an examination of the flux- 
dependence of the variability. Below we discuss the results and, 
where possible, relate them to the known behaviour of GBHs (for 
which we use a fiducial mass of Mbh = 10 Mq). 



7.1 Power spectral shape in relation to other sources 

The power spectrum of the variations is well described by a bend- 
ing power law model that bends from an index of ai ow ~ —1.1 
to ahigh ~ —2.3 around t^bend ~ 2 x 10 -4 Hz. This is fairly 
consistent with the lM c Hardv et al] J2006h relation 

log(Tb) = 2.1 log(Me) - 0.98 logfX^) - 2.32, (6) 

where Tb is the break timescalq'f|in days (i^bcnd = 1.16 x 10~ 5 /Tb 
Hz), Mbh = M 6 10 6 M , and L Bo i = L44IO 44 erg s" 1 is the 



5 Here we follow the notation used bv lKashvap et alj j20ld) and use /3 to 
denote the power of the test, i.e. the probability of meeting the detection 
criterion assuming the effect is real (e.g. there is a QPO present). 



6 Note that M c Hardv et at] fe006l) denned T b as the break frequency es- 
timated by fitting power spectra with a sharply breaking power law model. 
This difference is of little consequence to the present discussion as we ob- 
tained very similar estimates for the characteristic timescale using smoothly 
bending and sharply broken models (sectionfj). 
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10" 11 erg s _1 
3.5 x 10 41 erg s _1 at an assumed 
10 43 (assum- 



bolometric luminosity. The average 2 — 10 keV X-ray flux dur- 
ing the 2009 observations was F2-10 ~ 1.2 x 
cm -2 which gives L2 -10 
distance of 15.7 Mpc jRussellll2002l) . and Lboi 
ing a Z/Boi/^2-10 ~ 27; Elvis et al] [T994). The bl ack hole mass 
is tho ught to be M BH ~ 1.7 ± 0.5 x 10 6 M^ jDennev et alj 
2009), therefore riiEdd ~ ^Boi/^Edd ~ 0.05. Together these 
give a predicted bend frequency of 8 x 10~ 5 Hz (Tb ~ 0.14 
day), within a factor ~ 3 of the observed value. In fact the es- 
timated 1 4) end from the 2009 da ta is closer to the value expected 
from the |M c Hardv et alj j2006h relation than that estimated by 
lM c Hardv et alj d2004|). when calculated usin g the revised black 
hole mass estimate dPennev et alj2009Tl2010h . 

The high frequency power spectrum above the bend rarely 
provides strong, unambiguous signatures of the accretion state 
of GBH systems, which tend to differ more at lower frequen- 
cies where features such as low frequency QPOs, flat-topped 
noise and broad Lorentz i ans o c cur (see e.g. Ivan der Klia 120061 : 



iMcClintock & Remillardl 120061 ; iBellonil l201tt ). lM e Hardv et alj 
1 2004) suggested NGC 405 1 may be analogous to a soft state (par- 
ticularly like that seen in Cygnus X-l) on the basis of the low fre- 
quency noise (a ntt -1 spectrum extending over two decades or 
more), the high frequency slope (ahigh < —2) and the value of 
z'bond (which tends to be higher in the soft state, making the ratio 
T b /M B H smaller). 

For the 2009 data we found a very similar same high fre- 
quency shape and therefore is consistent with this identification. 
Ignoring any m dependence and simply scaling frequencies in- 
versely with mass (fbcnd oc 1/A/bh), the bend frequency for NGC 
4051 correponds to ~ 34 Hz in a typical GBH. This is close to the 
i^bond ~ 23 Hz estima t ed fro m a soft state observation of Cygnus 
X-l by |M c Hardvetal.l ( l2004h . 

The state analogy is supported by radio observations. At 
the highest available resolution the source is resolved into a 
core and two oppo sitely directed sources separated from the core 
by ~ 0.4 arc-sec dUlvestad & Wilsoi3ll984 iKukula et all [19951 
iHo & Ulvestadll200ll ; lGiroletti & Panessall2009l) . This structure is 
similar to the core and hotspots commonly seen in Fanaroff-Riley 
II radi o galaxies, althoug h no radio jet has been detected in NGC 
405 1 . Ijones et all d2010h recently used radio monitoring observa- 
tions and found no evidence for strong variabilit y of the co r e flux . 
Pespite the presence of these radio sources, Ijones et 
showed that NGC 4051 falls below the 'fundamental plane' con- 
necting the radio and X-ray lumin osities to black hole mass in radio 
loud AGN and hard-state G BHs dMerloni et al]|2003l ; iFalcke et al .1 
120041 ; iKordingetalJ 12006). i.e. it is less radio-loud than would 
be expected for a hard-state or radio-loud system (see Fig. 14 of 
I Jones et al J2010h . This is consistent with the 'soft state' interpreta- 
tion of NGC 405 1 mentioned above. 

The XMM-Newton data reveal the power spectrum of NGC 
4051 up to frequencies as high as > 10~ 2 Hz (see Fig[2]and[4]l, 
a factor ~ 40 above the bend frequency z/bcnd- Assuming a black 
hole mass as above, the gravitational radius is r 9 = GM/c 2 w 
2.5 x 10 6 km, and the period of the innermost stable circular orbit 
(ISCO) should be ~ 100 — 8 00 s, depending on the dimension- 
less spin parameter, j (see e.g. Ivan der Klisll2006l. and references 
therein). The power spectrum of NGC 4051 above i^bcnd shows 
no indication of deviating from a power law up to the highest fre- 
quencies (but see section[73}, and in particular shows no signature 
of the ISCO period, which might be expected to produce an ex- 
cess of variability pow e r due to Popple r boosting of "hot spots" 
jRevnivtsev et al]|2000l ; Isvunvaevlll973l) . There is no sign of the 



additional power spectral break seen at very high fre quencies in 
the hard state of Cygnus X-l bv lRevnivtsev et all l feOOo fl. Indeed, 
assuming linear scaling of timescales with black hole mass, we de- 
tect power at frequencies corresponding to ~ 1.7 kHz for a GBH, 
far above the fastest variations yet dete cted from these sources 
dRevnivtsev et alj200d : lstroh maver 200ll) . 



7.2 Stationarity and the existence of "states" 

The location of the power spectrum be nd frequency is lower than 
that reported by M c Hardv et al. (2004) from their analysis of the 
2001 XMM-Newton data. The 90 per cent confidence intervals from 
these two analyses - [5, 12] x 10~ 4 mHz and [1.4, 3.3] x 10~ 4 Hz 
- do not overla p, although some of t his may be due to the inter- 
val reported by M c Ha rdv et alj |2004) being an underestimate (see 
Mueller & Madejski 2009). The difference between the two results 
may in part be due to differences in the data reduction and power 
spectrum estimation techniques; the re-analysis of the 2001 data 
presented here (section[6]l gives a lower ^bend estimate and shows 
the difference between the two observations to be only moderately 
statistically significant (compare also the Bayesian posterior dis- 
tributions of Fig[3](. If this is a real effect, it is different from the 
expected scaling of ^bend with accre tion rate in the central region 
implied by the liyTHardv et alj 12006) relation, since the X-ray lu- 
minosity in the month leading up to the 2009 observation (based on 
the RXTE monitoring) was ~ 20 per cent higher than in the month 
leading up to the 2001 observations, but the bend frequency was 
~ 60 per cent lower. 

The separation of ~ 8 years between the two observations cor- 
responds to a separation of only ~ 30 minutes for a 10 MqGBH; 
on these timescales GBH power spectra are approximately station- 
ary except during very rapid state transitions. By the same linear 
scaling each of the ~ 40 ks XMM-Newton observations is equiva- 
lent to ~ 0.2 s, and the 45 day span of the observations covers the 
equivalent of only ~ 20 s for a typical GBH. Unless NGC 4051 
is drastically more non-stationary that most GBHs we might there- 
fore expect all the XMM-Newton observations to be representative 
of almost exactly the same state and variability properties. 

Puring both the 2001 and 2009 observations the variability 
does, to a very good approximation, follow a single rms-flux re- 
lation over the full flux range (see Figure [5) and the shape of the 
power spectrum remains roughly the same at different flux levels 
(see Figure |4). The lack of variability during low flux intervals 
(e.g. revl739) and the extreme variability during high flux intervals 
(rev 1730) are opposite extremes of this relation. These apparently 
different behaviours are therefore not the result of distinct variabil- 
ity stated d espite their very d ifferent energy spectra (see Figure 
U. (See also lUttlevetai]|2003l .) The variations in energy spectral 
shape are also continuous, as discussed by Vaughan et al. (2010, in 
prep.). 



7 Above the break the power spectrum of NGC 405 1 is sufficiently steep, 
a high < — 2, that the integrated power in the power spectra of Lx and 
dLx/dt both converge at high frequencies. The rate of change of the lu- 
minosity may be limited by the radiative efficiency of accretion I Fabian 
Il97l . The power spectrum is therefore not required to become steeper at 
even higher frequencies. 

8 Not in the sense normally applied to GBHs where different states 
are often characterised by distinct power spectra l shapes and ampli - 
tudes as well as different e n ergy spectra. See [van d er Klisl ( 2006 ) ; 
iRemillard & McClintockj 12006); Mc Clintock & Remillardl fe006l) . 
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7.3 Absence of QPOs 

We found no strong evidence (using a a = 0.01 test) for QPOs 
(additional, narrow components) in the 2009 power spectrum, and 
we place limits of 7? < 2 per cent rms on the strength of a QPO 
in the > 3 x 10~ 3 Hz range, and < 5 — 10 per cent over the 
10~ 4 -3x 10~ 3 Hz range. The 10" 3 -0.1 Hz bandpass over which 
these data are sensitive to weak QPOs corresponds to > 170 Hz for 
a GBH, where HF QPOs are often found, and the upper limit of 
R ~ 2 per cent on th e QPO strength is com parable to the strength 
of typical HF QPOs l lRemillard etai]|2003l) . Such HF QPOs been 
detected in GBHs onl y in interme diate states (particularly the soft 
intermediate state; IB ellonil [20 1 Oh . but are not always present (at 
least not at detectable levels), and when detected may be strongly 
energy-dependent. The lack of detection in NGC 4051 therefore 
does not constitute strong evidence against an analogous state in 
this AGN. 

Despite not being formally significant in an a = 0.01 test, 
the best-fitting narrow Lorentzian QPO component has a frequency 
of vq = 4 x 10~ 3 Hz and rms R ~ 1 per cent. This frequency 
is very close to that found previously by Vaughan &Uttlev (2005) 
from fitting the power spectrum of the 2001 data. This may be a 
coincidence in the random sampling fluctuations around a smooth 
continuum power spectrum. But it is around this frequency that the 
signal-to-noise in the data is highest, and so model fitting will be 
most sensitive to small (~ 20 per cent) deviations from the smooth 
power law model. The high frequency power law may be only an 
approximation to a more complex continuum power spectrum, the 
structure of which is on the edge of detectability in these data. 

It is well known that HF QPOs in GBHs can be strongly en- 
ergy dependent, and frequently stronger at higher energies, and in 
fact we did find the higher energy band (2-10 keV) gave slightly 
stronger evidence for a possible QPO. But this could also indicate 
that this energy band contains a stronger contribution from a sec- 
ond, variable emission component with a different power spectrum 
from that which dominates the softer energy bands. 

7.4 A quasi-constant emission component 

The flux offsets in the energy-resolved rms-fiux relations, C(E), 
reveals the spectrum of a quasi-constant component (QC) of the 
emission. A simple model comprising a single power law (r = 
3.06 ± 0.13) modified by Galactic absorption felvis et alii 19891) 
gave a very good fit, with \ 2 = 10.65 for 18 dof. A Comptonising 
plasma model (compTT) with optical depth r = 1 and seed photon 
energy ~ 50 eV, which has an approximately power law shape, also 
gave a good fit (% 2 = 9.45 for 17 dof). The simplest interpretation 
is that this represents an emission component that does not vary at 
all during the 2009 observations; emission that remains if the highly 
variable emission component(s) drop to zero flux. We note that the 
Chandra data taken during an extended low flux period in 2001 
and the XMM-Newton observation taken during a low flux period in 
2002 were previously interpreted in terms of an unresolved (<; 100 
pc), non-varying soft spect ral component and a s trongly variable, 
harder spectral component dUttlev et alj2003l , l2004l) . However, the 
QC emission is not required to be absolutely constant over the ob- 
servations, it is only required to have a far lower (absolute) ampli- 
tude (at > 1 mHz) and be weakly correlated (or uncorrelated) with 
the highly variable component that dominates brighter flux spectra 
(e.g. revl730) and drives the rms-fiux relation. 

The spectral shape and absolute flux level of this component 
are similar to the lowest flux spectrum taken during 2009 (revl739) 



below ~ 2 keV, as shown in Figure [6] This strongly suggests that 
most of the soft X-ray continuum emission observed during low 
flux periods is due to the soft QC component. This may be the tail 
of the spectrum of inverse-Compton scattering of optical-UV con- 
tinuum photons in shock-heated gas where the ionised outflow col- 
lides with the inter-stellar medi um of the h ost galaxy, as recently 
predicted to occur in AGN bv lKind ([20 10). Indeed, the velocity- 
ionisation structure of the outflow as resolved by the R GS data pro- 
vide further evidence to support this scenario JPounds & Va ughan 

[ml- 
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APPENDIX A: CONSTRUCTION OF 'FLUXED' SPECTRA 

X-ray spectra are conventionally given in terms of count rates per 
detector channel. These "raw" spectra are of limited use for data 
visualisation because of the distorting effect of (i) the energy- 
dependent efficiency of the telescope and detector (described by the 
effective area function), and (ii) blurring by the line spread function 
of the instrument (modelled by the redistribution matrix). The true 
source spectrum S(E) (in units of photon s _1 cm~ 2 keV -1 ) is re- 
lated to the observed count rate in channel i by a Fredholm integral 
equation: 

C(i) = J R(i,E)A(E)S(E)dE (Al) 

where A(E) is the instrument effective area (cm 2 ) and R(i, E) 
gives the probability that a photon of energy E will be recorded in 
channel i (i.e. the redistribution function). Here the observed spec- 
trum C(i) is assumed to be background-subtracted, and the instru- 
mental response defined by R(i, E)A(E) is assumed to be linear, 
that is, independent of the value of S(E) which is a good approx- 
imation so long as pile-up is negligible. This integral equation can 
be approximated using a discrete formulation in which the energy 
range is divided into iV bins Ej : 

N 

c, \J /,',.. t,.s; ; (A2) 

2=1 

where Rij is the average and Aj and Sj are the definite integrals 
of their continuous counterparts, over the energy range of each bin 

Ej E J+1 . 

In general, methods that attempt to invert the above equations 
to recover S(E) from an observed spectrum d are unreliable and 
one must turn to the 'forward fitting' method to build a model of 
the spectrum (see Arnaud 1996). However, it is often useful to view 
even an approximate 'fluxed' spectrum, largely free from the dis- 
torting effects of the instrument efficiency, if only to giv e a bet- 
ter vis ual representation of the underlying source spectrum. iNowald 
(2005) discussed a simple method to accomplish this in which the 
observed count spectrum is normalised by an effective area curve 
that has been blurred using the redistribution function. 

JV 

Si = C i /^2R ij Aj (A3) 
i=i 

This method matches the blurring of the effective area to the blur- 
ring of the observed spectrum to give an estimate of the spectrum in 
true flux units (e.g. photon s _1 cm~ 2 ) that is, to a large extent, free 
from the distorting effects of the energy-dependent effective area. 
This process is essentially equivalent to calibrating the sensitivity 
using an observation of an ideal standard star (with a perfectly flat 
spectrum|f|. The resulting spectrum can easily be converted to con- 
ventional flux density units (e.g. erg s _1 cm" ~ 2 keV" 1 ). 

It is crucial that the effective area function Aj is blurred by the 
redistribution matrix. The reason is that the effective area function 
may contain abrupt changes in efficiency (e.g. near photoelectric 

9 This can be done trivially within XSPEC by defining a constant model 
spectrum M(E) = 1.0 (e.g. a power law with index and normalisa- 
tion 1) and plotting the "unfolded" spectrum. For the special case of a con- 
stant model this is equivalent to equation |A3l The 'unfolded' plot shows the 
model in flux units (= 1) multiplied by the ratio of the observed spectrum to 
the 'folded' model, i.e. M(E;) x d/ J2j Rij A o M j = Gil Ej Rij A j- 
the SAS task ef luxer can also perform this transformation. 



edges in the detector-mirror system, such as the K-edges of O and 
Si, or the M-edge of Au) that can be far sharper than the resolv- 
ing power of the detector, which means that normalising the ob- 
served (blurred) spectrum by the exact (not blurred) effective area 
introduces very strong spurious features into the spectrum near the 
edges. Blurring the effective area curve matches the resolution to 
that of the data and suppresses this effect. The result is a spec- 
trum with the energy-dependence of the detector/mirror efficiency 
removed (to a reasonable approximation), but which remains at the 
detector spectral resolution. We emphasise that this is not, even ap- 
proximately, a deconvolved spectrum. It has the advantage over the 
popular 'unfolding' method (in XSPEC) of being independent of 
the spectral model, and therefore providing a more objective visu- 
alisation 'fluxed' spectrum. 



APPENDIX B: QPO DETECTION AND UPPER LIMITS 

This appendix presents first an outline of the method used to search 
for QPOs in the power spectral data, and second how upper limits 
were derived based on the detection method. 



Bl Calibrating the QPO detection procedure 

An a-threshold detection procedure uses a test statistic T to quan- 
tify evidence for a detection, and is calibrated such that a detection 
occurs when T(x) > T cr i t for some threshold value T cr j t cho- 
sen to have a small probability a of occurring by chance when 
Hq is true. Specifically, the detection threshold is chosen such that 
Pr(T > T cr it | Ho) ^ ot, where a is the chance of a 'false de- 
tection' (a 'type I error' when Ho is true). Commonly used values 
are a = 0.05 and 0.01. But in order to calculate T cr it we need to 
know the reference distribution of the test statistic on the assump- 
tion that the null hypothesis is true, p(T \ Ho). In general this 
may be estimated by calculating the distribution of the test statistic 
from a large number of Monte Carlo simulations of data generated 
assuming Ho and distributed as expected for realistic data. 

The above would be valid if the null hypothesis was sim- 
ple, that is, did not contain any unknown (free) parameters. In 
the present case the null hypothesis does involve parameters es- 
timated using the observed data, 6c- We account for the uncer- 
tainties in their values by using random parameters drawn from 
the posterior distribution to generate the simulated data, using the 
MCMC discussed in section[3] This gives the distribution of t he test 
statistic known as the posterior predictive distribution (see Rubin 
ll984l:lMenelll994l : lGelman et al.ll996ll2004|Protassov et alj2002T 
IVaughanll2010l) . 

P (T | Ho, x obs ) = J P (T | Ho, 6 c )p{ecW hs )d0c (Bl) 

In order to determine T cr i t this distribution was mapped out using 
Monte Carlo simulations as follows: 

• repeat for k = 1,2, N simulated datasets: 

(i) Draw a set of parameter values from the posterior 9 C 

(ii) simulate data from Ho with parameters 6 C : x fc 

(iii) fit# togetx 2 (x' ; ,Jfo) 

(iv) fittfitogetx 2 (x fc ,#i) 

(v) compute T h 

• use the distribution of T k to define T crit for which Pr(T > 

T crit I Ho) < a 
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We used TV = 3, 000 simulations to calibrate T cr it and found 
the critical value for an a = 0.01 test was T cr i t = 12.28. (This 
was calculated using a Q — 10 QPO, but using Q — 5 gave a 
very similar result.) for the data and models discussed in this pa- 
per. In order to ensure the global minimum was found in each case 
(this can be tricky for Hi which may produce multiple minima 
in y 2 ) we used the robust automated fitting scheme discussed by 
lHurkettetai] ( l2008l. section 3.2.2) 

B2 Upper limit procedure 

We wish to determine the upper limit U on the QPO strength R, 
which is the smallest value of R such the probability of detecting 
the QPO (assu ming it is present, i.e. Hi) is ^ /3 m in for some spe- 
cific Pmin- See lKashvap et alj d2010h for a very clear introduction 
to the meaning and calculation of upper limits in astronomy. 

In the present case we are searching for a QPO of unknown 
location (frequency) but the upper limit may vary as a function of 
this location - because the "background" continuum level and the 
uncertainty on the spectrum vary as a function of frequency - so 
we must calculate the limit for a set of different QPO locations vq. 
The probability of a detection, ft, which is 1— (probability of 'type 
II error'), is: 

p(v Q , R) = Pr(T > T crit | Hi,v Q ,R) (B2) 
and the upper limit U(uq) is the smallest value of R which satisfies 

P(y Q ,R) Z/3 mln (B3) 

for some chosen minimum power /3 m m- 

Given the T cr i t defined above we can compute the upper limit 
using Monte Carlo simulations of realistic data under Hi. These 
can be used to map out the distribution p(T | Hi, vq, R) for dif- 
ferent vq and R values, and hence find U(vq) that satisfies the 
above equations. In these terms the upper limit is the strength of the 
weakest QPO that would have a reasonable probability ft ^ fimirx, 
for some specific fi m in such as 0.5 or 0.9. 

This method will, in general, be sensitive to the parameters of 
the continuum model, 9c, which of course we do not know exactly. 
But, as before, we can average over the posterior distribution of 
these as derived from the data p(9c | x obs ). 

P (T | Hi,v Q ,R) = J p(T\ Hi, v Q , R, 6 c )p(e c j x° bs )d0 c . 

(B4) 

Again, this can be done be tacking values from the posterior using 
the MCMC discussed in section[3] 

The calculation of ft for a given pair of QPO parameters (vq 
and R) was performed using a Monte Carlo scheme as follows. 

• loop over k = 1, 2, TV simulated datasets 

(i) Draw a set of parameter values from the posterior 6q 

(ii) simulate data from Hi with continuum parameters 8 C and 
QPO parameters (vq,R): x fc 

(iii) fitffotogetx^x*,^) 

(iv) fittf 1 toget X 2 (x fc ,tfi) 

(v) compute Tk 

• use distribution of T k to find P(vq,R) = Pr(T > T cri t | 
Hi,v Q ,R) 

This process is repeated over a grid of QPO parameters, 
vq and R, and at each frequency vq the smallest R for which 



P(vq, R) Pmin is the upper limit U (vq). For the present analy- 
sis we used 30 logarithmically spaced values of QPO location vq, 
40 logarithmically spaced values of QPO strength R, and at each 
pair generated TV = 200 simulations. This is a sufficient number 
of simulations to accurately define contours of e.g. /3 = 0.9 on the 
vq-R plane. 

The detection procedure was calibrated using simulated data 
generated assuming Ho. But to evaluate /3, and hence the upper 
limit, the data were simulated assuming Hi. The former is con- 
cerned with the chance of spurious detection assuming no QPO 
present, the latter is concerned with the chance of detection when a 
QPO is present. 

This paper has been typeset from a TgX/ ETpX file prepared by the 
author. 



